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Abstract. In a first course of general relativity it is usually quite difficult for students to grasp 
the concept of a geodesic. It is supposed to be straight (auto-parallel) and yet it 'looks' curved. 
In these situations it is very useful to have some explicit examples available which show the 
different behaviour of geodesies. In this paper we present the GeodesicViewer, an interactive 
tool for studying the behaviour of geodesies in many different space-times. The geodesies can 
be represented in several ways, depending on the space-time in question. The use of a local 
reference frame and 'Cartesian-like' coordinates helps the students to develop some intuition 
in various situations. We present the various features of the GeodesicViewer in the form of 
readily formulated exercises for the students. 
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1. Introduction 

The intrinsic curvature of a space-time in general relativity is a concept that contradicts our 
every day experience of space and time. The necessary mathematical set of tools is difficult 
to learn and fairly abstract. To get some impression what a curved space-time means, is to 
study the behaviour of light rays and particles in free motion. In the geometric optics limit 
and for particles whose mass has no back-reaction on the curvature of space-time, light rays 
and particles in free motion can be represented by null and time-like geodesies, respectively. 

For a first glimpse on how null and time-like geodesies behave, it may be sufficient to 
use off-the-shelf/standard software like for example Maple, Mathematica, or Octave. All of 
them could integrate the geodesic equation and show the geodesic as 2d- or 3d-plot. However, 
to explore the behaviour of geodesies, an interactive tool is indispensable. 

In this article we present the GeodesicViewer an interactive tool to thoroughly 
examine the behaviour of light-like and time-like geodesies in a space-time whose metric 
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Figure 1. Screenshot of the GeodesicViewer's user interface. 



is provided analytically. The database of metrics is taken from the Motion4D library ID. The 
metrics with the corresponding Christoffel symbols and local tetrads are detailed in Ref. l3\. 
The graphical user interface, see Fig. [T] is written using the object-oriented, cross-platform 
application framework Qt[41. The graphical 2D and 3D output is realized by means of the 
Open Graphics Library (OpenGL)15i. For the numerical integration of the geodesies, we use 
a standard fourth-order Runge-Kutta method and the integrators of the GNU Scientific Library 
(GSL)|6|. 

The structure of the paper is as follows. In Sec. |2]we give a short description of the 
GeodesicViewer. Sec. [3] discusses the standard situations in the Schwarzschild space-time. 
Periodic orbits of time-like geodesies in black hole space-times are worked out in Sec. |4] 
Secs.|5]and|6]deal with the more exotic space-times like the Morris -Thorne wormhole and the 
extreme Reissner-Nordstr0m diblack hole. 

The GeodesicViewer is freely available for Linux and Windows. The source code and 
several examples can be downloaded from www.vis.uni-stuttgart.de/relativity. 



2. GeodesicViewer 



The two main outputs of the GeodesicViewer are 3D and 2D representations of the geodesic 
data. In the standard 3D representation, a geodesic is depicted by means of pseudo-Cartesian 
coordinates, where the inherent coordinates are transformed into Cartesian coordinates as 
usual. If an embedding diagram is defined for a specific hypersurface of the space-time, 
the geodesic can also be represented in this form. (In principle, the user can implement 
any representation he wants.) The standard 2D representation of a geodesic is also given 
in pseudo-Cartesian coordinates, where now only a specific hypersurface is used. Another 
representation shows coordinate or velocity relations. For example, the radial coordinate 
could be plotted against the affine parameter. A third representation follows from the Euler- 
Lagrangian formalism, where an effective potential can be defined. This representation is 
particularly helpful to find bound orbits. 

The numerous features of the GeodesicViewer are described by means of examples which 
are formulated as exercises. The mathematical explanations are directed to the teacher. The 
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exercises can be made either from scratch or a configure file can be prepared in advance, so 
that the student only has to change a few parameters. Each result description is accompanied 
by a configure file which holds the final result for reproduction in the GeodesicViewer . To 
familiarize oneself with the graphical user interface of the GeodesicViewer, there are several 
online tutorials. These can be worked through either alone or with guidance from the teacher. 



3. Basic examples in the Schwarzschild space-time 

The prime example of general relativity is the Schwarzschild metric which we will give here 
in isotropic coordinates — {t,x,y,z). The line element ds^ — g^^^dx^dx'^ reads 

ds" = (^i-^^c^^f^^ [dx^+dz+d^] , (1) 

where = x^ +y^+z^, Ps = GM/ (2c^) is the Schwarzschild radius, G is Newton's constant, 
M is the mass of the black hole, and c is the speed of light. The transformation between 
the usual Schwarzschild radial coordinate r and the isotropic radial coordinate p is given 
by /- = p (1 +ps/pf. If M = 0, Eq. Q simplifies to the Minkowski metric in Cartesian 
coordinates. 

Because of the spherical symmetry of the Schwarzschild space-time, we can restrict 
geodesies to the xy-plane. Then, the geodesic equations read 

o = f+2r;,fi+2r;,,f3), (2) 

= X + nf + rj,i2 + 2T'^.xy + r^^f , (3) 
0^y + T^,f + TlJ^+2r%xy + Ty^, (4) 
with the Christoffel symbols 

_ 2c^p3p,(p-p,)x _ 2c^p^p,{p-p,)y 

iP+Ps) (P+P.v) 

, _ 2p,x , _ 2p,y 

p3[l-p2/p2]' ''y p3[l_p2/p2]' 

=v =-p' ^_2p, X 

" ■^>' p3 i+p,/p' ^ ' 

r^' — — r^' = r^' — — (q\ 

»■ p3 l+p,/p' ^ ' 

Here, a dot represents the derivative with respect to the affine parameter X, hence i = dt/dX. 

To integrate the geodesic equations, we need not only an initial position but also an 
initial direction. For this purpose, we first introduce the local reference frame {e(,) },=,.,;,,■ ^ of 
an observer which defines a local Minkowskian system. The four base vectors e(,-) — e|^j(9^ 

have to fulfill the orthonormality condition ^^vej^^e^^j ^^'^ ^(')(;) = ^^^^§("1' 1' 

Here, the most convenient choice for the local reference frame is the one which is adapted to 
the coordinates and the symmetries of the metric. 



\+Ps/pd, ^ / , p 



--l-p.,/pc' e,,= ^l + ^j\, (9) 

+ '«^0+p) ^^^^ 
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Figure 2. Initial direction k = i''' e^,) with respect to the local reference frame of an observer. 

Now, an initial direction, k = k'^'^^e^.^dfi — k^d^, of a null or a time-like geodesic can be defined 
with respect to this local reference frame as shown in Fig. [2] 

For a null geodesic we have k — ejoj + sin jcos^e(i) + sinxsin^e(2) + cosx^{3)- The 
initial direction of a time-like geodesic, on the other hand, equals the initial four-velocity 
k = u = yccjo) + 7j3c (sin J cos 1^6(1) + sin j sin ^6(2) +cos je^)), where j3 = v/c is the initial 
velocity v scaled by the speed of light, and y = l/i/l — j3^. Because we restrict to geodesies 
in the xy-plane, we must set X = 90°. The geodesic equations can now be solved using the 
initial position -^^ and the initial direction q = k^. 

The numerical integration of the geodesies are accomplished by either a standard fourth- 
order Runge-Kutta method or by the integrators of the GNU Scientific Library. After each 
integration step, the geodesic is tested if it still fulfills the constraint equation gfivx^x"^ — K"c^ 
with K — QfoT light-like and fc = — 1 for time-like geodesies. If this constraint is not fulfilled 
within a certain accuracy, the integration stops. 

For the following examples, we set G = c = M = 1. In physical units, this means that 
G/c^ = 1 Isj ^ which equals one light second per mass unit ^ « 2.03 x IO^Mq. Hence, 
M = 1 represents an object of 2.03 x 10^ solar masses and distances are measured in light 
seconds. 

3.7. Deflection of light 

In the neighbourhood of a black hole, light rays are no longer straight lines but they are 
deflected due to the curvature of space- time. The size of the deflection depends on how close 
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the light ray passes the black hole. 

Exercise: Start the GeodesicViewer and select the Schwarzschild metric in 
isotropic coordinates (Schwarzschildlsotropic) in the "Metric/Integrator/Cons- 
tants" window. Set the position of an observer in the "Local Tetrad" window 
to (f = 0,x = 6,y = 6,z = 0). With respect to his local reference frame, the 
observer starts light rays in the negative y-direction, k = e(o) ~ ^(2)- For that, 
set the direction angle B, — 180° within the "Geodesic" window. Modify the 
initial position y and describe the behaviour of the light rays. What is the closest 
approach where the light ray can still escape from the black hole? 

Result: Figure |3] shows light rays for several initial positions y,. The closer the 
light ray passes the black hole the stronger the rays are deflected. If the light ray 
crosses the photon orbit, r^ ~ |rj (Pp„ = (2 + -\/3)p.s), the light ray cannot escape 
from the black hole. 




Figure 3. Light rays starting from a: = 6, y, = {6,5.5,5,4.5,4, 3.5, 2} in the negative y-direction 
are deflected by the curved space-time close to the Schwarzschild black hole (black disk). The 
dashed circle represents the photon orbit which represents the limit of the closest approach for 
light rays. 



Configure file: defl_of_light. 



A detailed explanation of the configuration file can be found in Appendix A 



It could also be interesting to find the initial position y where the light ray approaches 
the photon orbit asymptotically. Here, we have y k, 3.9508265. For that, however, it would 
be better to change the numerical integration from the standard Runge-Kutta fourth order 
integrator without step size control to the Runge-Kutta-Fehlberg integrator with step size 
control and absolute error tolerance of e^bs — 10^'^. Then, the constraint equation for the 
light ray is fulfilled until the light ray reaches the horizon. 



3.2. Einstein ring 

In the previous example, we have seen that for some specific initial positions y, the light ray 
orbits the black hole and returns to the observer Due to the symmetry of the Schwarzschild 
space-time, this happens also for light rays which lie in a plane that is rotated around the 
connecting axis between the observer and the black hole. As a result, the observer will see 



Studying null and time-like geodesies in the classroom 



6 



himself as an Einstein ring around the black hole. However, there is not only one but an 
arbitrary number of Einstein rings. 

One method to find the geodesies that return to the observer would be to solve the 
geodesic equations analytically, see e.g. Miiller |7|. But this is extremely laborious even for 
the more adequate spherical coordinates. Another possibility, known from solving ordinary 
differential equations, is the shooting method which we can be easily simulated using the 
GeodesicViewer and varying the initial direction ^. 

Exercise: Start the GeodesicViewer and select the Schwarzschild metric in 
isotropic coordinates (Schwarzschildlsotropic) in the "Metric/Integrator/Cons- 
tants" window. Set the position of an observer in the "Local Tetrad" window to 
(f = 0,x ~ 6,y = 6,z = 0). Find the initial directions t, of the light rays that orbit 
the black hole once, twice, or thrice, cf. Fig. B] before returning to the observer 




Figure 4. Light rays starting from the observer at (x = 6, v = 6) return to him after orbiting 
the black hole once (solid line) or twice (dashed line). Hence, the observer will see himself as 
an Einstein ring of first or second order. 



Result: The initial directions ^ with respect to the observer's local reference 
frame, cf. Fig. |2] read: E,i = 195.668165°, ^2 = 195.96490458°, and ^3 = 
195.965451625°. 

Configure file: einstein_ring_l,einstein_ring_2,einstein_ring_3. 

Exercise: The situation is as before. But now, the observer sends out a flash 
of light and measures the time until the Einstein rings appear: X\ ~ 40.1480, 
T2 ~ 69.2052, T3 ~ 98.2203. What is the relation between proper time and 
coordinate time? Switch to the 2D view and use the coordinate time 'f ' as abscissa 
and 'x' as ordinate. 

Result: From Fig. |5] we can read the light travel times f, for the light rays with 
initial directions ^. The relation between the coordinate time f, and the proper 
time T; for a fixed observer position is constant. 

Configure file: lightJraveLtime. 

This example could be also formulated for given observer positions p and corresponding light 
travel times T(p) as measured by the observer To prepare this exercise, the coordinate times t 
for one orbit have to be determined in advance for several radial positions p. Then, by means 



of the relation x/t = yj\ -r^jr = w 1 -4ps/ p (1 +p.s/p) 



^2 



which shall be found by the 
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Figure 5. Light travel times for the light rays of Fig. [4] in 2D coordinate representation: 
(1 «45.1758,r2 - 77-8718, f3 110.5205. 

Students, the proper times T can be determined. Probably, this relation could be easier found 
in the standard spherical Schwarzschild coordinates. 

3.3. Shapiro time delay 

In Ref. [8], Rindler discusses the Shapiro time delay by means of a situation similar to the 
one shown in Fig. |6] A light ray emitted at x„,„,„ = —X passes a spherical mass M, which is 
located at (x = 0,y = 0), ?Ay = R and reaches the observer at x„fe = X. From the isotropic 
metric, Eq. together with ds = and p = \/X'^ +R^, he deduces the total coordinate time 

Af «2X + r,ln-^, (11) 

where the logarithmic term is the Shapiro time delay. 

Exercise: Start the GeodesicViewer and select the Schwarzschild metric in 
isotropic coordinates. Find a geodesic that connects the observer, located at 
x„bs — 200, with a light source at position x„„„ = —200. Compare the total light 
travel time with the approximation by Rindler 




Figure 6. Shapiro time delay. The solid line represents the real light ray, whereas the dashed 
line is the rectilinear approximation. Impact parameter: y = R. Position of observer and source: 

^obs — ^source — 



Result: The light ray that connects the observer and the light source in Fig. |6] 



has initial angle ^ « 173.8578332°. Rindler's approximation, Eq. (jllb, yields 
At « 412.0797, whereas from the GeodesicViewer we obtain At « 414.579 

Configure file: shapiro. 
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4. Periodic orbits of time-like geodesies in black hole space-times 

The most simple periodic orbit of a time-like geodesic in the Schwarzschild space-time is a 
circular orbit. The relation between the radius r of the orbit and the corresponding velocity j3 
can be easily found within the original, spherically symmetric Schwarzschild coordinates. In 
isotropic coordinates, this relation reads 



i3-' = 



\ ^ [4p.s V' ' P 



(12) 



which is valid only down to the last stable orbit, p > p,„ — (5 + 2\/6)Pi. However, beside 
these circular orbits, there are also more complex periodic orbits. Levin and Perez-Giz |9l 
give a whole taxonomy of periodic orbits. See [Appendix B| on how to reproduce the orbits 
shown in their paper. 

In the following exercise, we will consider periodic orbits for Kerr black holes. In Boyer- 
Lindquist coordinates, see e.g. Bardeen et al. ifTOl . the Kerr metric reads 

ds^ = -{\-^] c'-dt^ ^— cdtd(p + -dr^+T.d-&^- (13) 



E 



ryfl^r sin^ i> 



9 9 AjU r bill \ . 9 „ , 9 

r^ + a^ + - jsm^Mq)^, (14) 

' fl^cos^ jj, A = — r^r + a^, and rj — IGM/c^. M is the mass and a is the 
angular momentum per unit mass of the black hole. As the observer's local reference frame 
we use the locally non-rotating tetrad 



A / K . \ /A 



" V EA y^"' '^'""fj' " V E*^" ^^^^ 

1 1 
^ ' y/Y. V A smz> ^ 

where CO = VsarjA and A = (r^ + a^) E + r^a^r sin^ t?. 

Exercise: Select the Kerr metric in Boyer-Lindquist coordinates (KerrBL) and set 
the mass and angular momentum parameters to M = \, a = 1/2. The observer 
is located at (r — 6,1!} — 7t/2,(p — 0). 'With respect to his locally non-rotating 
reference frame, cf. 'Natural local tetrad' in the 'Local Tetrad' window, he starts 
a future directed time-like geodesic with initial velocity /3 = 1/2 and direction 
1 = 90°, ^ = 67.649°, see Fig.^ By varying the initial direction t,, find periodic 
orbits of higher order Instead of varying the initial direction, the initial velocity 
could also be changed. 

Result: Periodic orbits of higher order follow from initial directions ^2 ~ 82.95°, 
« 72.6°, ^4 ~ 70.58°, ^5 « 75.1°, etc. For the fixed direction ^1 = 67.65°, we 
can also change the initial velocity: ^2 ~ 0.545, jSs ~ 0.52006, jSs w 0.528. 

Configure file: kerr_periodic. 



5. Wormhole space-time 

The most simple non-trivial space-time is that of a Morris-Thorne ifTTl wormhole whose 
metric, given in spherical coordinates, reads 

ds^ = -c^dt^ + dl^ + {bl + f-) {d^^ + sin^^ d(p^) . (17) 
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Figure 7. Periodic orbits in tlie Kerr space-time witli parameters M = 1 and a = 1/2. The 
observer is located at r = 6, (p = and stalls a time-like geodesic with /3 = 1/2 and % = 0°, 
£,1 = 67.65° (left), or ^2 ~ 82.95° (right). 



Here, bo is the throat radius and I is the proper radial coordinate. The local reference frame 
of a static observer is given by 



1 



i(r) 



1 



-/2 



1 



- /2 sin t> 



(18) 



To receive an impression of how the inner geometry of this wormhole space-time looks 
like, we take advantage of the spherical symmetry and embed the (f — const, = n/2) 
hypersurface into the three-dimensional Euclidean space, cf. for example Ref. lfT2l . The 
resulting embedding function reads 



z{r) = ib^oln 




1 



(19) 



where the radial coordinates r and / are related via — bl + fi. Figurejsjshows an embedding 
diagram for a wormhole with throat size bo = I and an observer at (/ = 10, = ;r/2, (p =0). 



Studying null and time-like geodesies in the classroom 



10 



i = 14° 




Figure 8. Embedding diagram of a Morris-Tliome wormhole space-time with foo = 1 and some 
exemplary null geodesies (black lines) with initial direction k = e^,) — cos^e^;) +cosi^e(^) 
starting at the observer's position (/ = 10, i? = ll/2, ip = 0). 



Exercise: Select the Morris-Thorne wormhole and set the throat size parameter 
bo = 1. The observer is located at {I — — 7r/2, (j? = 0) and uses the inward- 
oriented tetrad e(o) = ^(t), ^(i) = ^^(l)> ^(2) = C(ip)> ^(3) = local reference 
frame. A light ray with initial direction % — 90°,(^o — 0° traverses the wormhole 
and hits the point P = (l — —10, (p — 0). Now, find the initial direction t,i where 
the new light ray hits the point P again but now with (p — 2%. Explain what an 
observer would see between and 

Result: A light ray with initial angle ^ « 5.64336° travels around the wormhole 
throat before reaching the point P. Between (^o and the observer would see the 
whole lower universe. Here, an interactive tool is indispensable to comprehend 
this fact. 

Configure file: morristhome. 

Detailed discussions for the first-person visualizations of the Morris-Thorne wormhole can be 
found in MuUer lUl or Ruder lfT?l . 

6. Multi black hole solution 

Interesting general relativistic situations occur when several black holes are combined. These 
configurations will be inherently dynamical due to the mutual gravitational attraction between 
the black holes. However, it is possible to construct static situations if one allows the black 
holes to carry an electric charge. Then the gravitational attraction can be compensated 
exactly by the electric repulsion. These solutions have been discovered by Majumdar ifTSll 
and Papapetrou |fT6l . They can be described as a collection of any number of extreme 
Reissner-Nordstr0m black holes. We consider here the case of only two such black holes, 
see Chandrasekhar |,17] . 

In Cartesian coordinates, the extreme Reissner-Nordstr0m metric reads 



where f/ = 1 +Mi /ri +M2/r2, r\ = ^/x^ +y^ + {z- 1)^, and = ^/x^ +y^ + {z+ 1)^. Here, 
we use geometric units with G = c = 1. The two black holes are located atx = y = 0,z = ±l. 




(20) 
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As local reference frame, we use 

e(,) = t/a, e(,) = e(,.) = l^y, e(,) = ^a. (21) 

For simplicity, we set M\ — M2 — 1 . 




Figure 9. Time-like periodic orbits in the extreme Reissner-Nordstr0m metric with the two 
black holes located at z = ±1. Left: The geodesic starts at {x = y = O.z = 3) with four- 
velocity u = 7(e(,) + /3e(;r)), where /3 0.591943 and 7= Xj^JX-fP-. Right: The geodesic 
starts at (y = z = 0,-V = 3) with four-velocity u = /(ej,) + j3 cos^ej,) + j3 sini^ej,,)) , where 
j3 = 0.2 and I, « 69.26. 



Exercise: Select the extreme Reissner-Nordstr0m metric (ExtremeReissnerNord- 
stromDihole) and set the mass parameters of the two black holes to Mi = M2 = 1. 
Find time-like periodic orbits in the xy- and xz-planes. 

Result: Some exemplary time-like orbits are the following: 1) {x — Q,y ~Q,z = 

3), u = r (e(,) + j3e(,)), j3 « 0.896899; 2) (x - 0,y = 0,z = 3), u - 7(e(,) + j3e(,)), 
J3 « 0.591943; 3) {x = 0,y = 0,z = 10), u = 7(e(,) + i3e(.,)), j3 w 0.451; 
A light-like orbit is given by {x = 0,y = 0,z~ 1.720218), k = e(,) +e.(x)- 

Configure file: extrRNdiholel, extrRNdihole2, extrRNdihole3, extrRNdihole4, 
extrRNdihole5. 



7. Summary and outlook 

The GeodesicViewer is a valuable tool to study the behaviour of light-like and time-like 
geodesies in space-times whose metrics are provided analytically. While geodesies could 
also be visualized using standard software, a clear understanding of their behaviour can only 
be achieved by interactively manipulating the corresponding parameters. 

In this paper, we have given a small number of possible applications formulated as 
exercises for students in the classroom. This can only hint at various other possibilities. 
The strength of GeodesicViewer lies in its visualization capabilities and to a large part in 
the available database of exact solutions. This allows the user to study the behaviour of 
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geodesies in many different circumstances. The possibility to choose the curve parameters in 
an interactive way means that one can make numerous parameter studies depending on the 
problem under consideration. 

So far, the GeodesicViewer can only handle one geodesic at once. This limitation will be 
overcome in a future version. We also plan to make the GeodesicViewer scriptable to realize 
more complex demonstrations. 



Appendix A. Configuration file example for the GeodesicViewer 

As an example, we give here the configuration file for the 'Deflection of light' exercise of 
Sec. 3.1 as written by the GeodesicViewer. Although this file is a plain text file, we do not 
recommend to modify it by hand. Besides, note that not all of these parameters are crucial for 
this exercise. 



METRIC 
PARAM 
INIT_POS 
INIT_DIR 
INIT_ANGLE_VEL 
TIME_DIR 
AXES_ORIENT 
GEOD_SOLVER_TYPE 
GEODESIC_TYPE 
STEPSIZE_CTRL 
STEPSIZE 
STEPSIZE_MAX 
EPSILONS 
CONSTR_EPSILON 
MAX_NUM_POINTS 
TETRAD_TYPE 
BASE_0 
BASE_1 
BASE_2 
BASE_3 
BOOST 

SPEED_OF_LIGHT 

GRAV_CONSTANT 

DIELECTRIC_PERM 



Schwarzschildlsotropic 
mass 1.000000000000 



6.00000000 
0.00000000 
90.00000000 



0.00000000 
-1.00000000 
180.00000000 

1 


4 

lightlike 
1 

l.OOOOOOOOe-02 
l.OOOOOOOOe+00 

l.OOOOOOOOe-12 . OOOOOOOOe+00 

l.OOOOOOOOe-06 

3000 



1.00000000 
0.00000000 
0.00000000 
0.00000000 
0.00000000 

1.000000 
l.OOOOOOe+00 
l.OOOOOOe+00 



6.00000000 
0.00000000 
0.99000000 



0.00000000 






00000000 





00000000 





00000000 


1 


00000000 





00000000 





00000000 





00000000 


1 


00000000 





00000000 





00000000 





00000000 


1 


00000000 


90 


00000000 





00000000 







Most of the parameters are self-explanatory. The time-direction (TIME_DIR) can be either 
plus or minus one depending on the geodesic being future- or past-directed. The method 
for solving the geodesic equation (GEOD_SOLVER_TYPE) is encoded by just a number, cf. 
m4dMotionList .h of the Motion4D library. The absolute and relative error tolerances are 
given by the EPSILONS values. The numerical integration of the geodesic stops when the 
constraint equation with CONSTR_EPSILDN is no longer fulfilled. For the local reference 
frame in this exercise, we use the default natural local tetrad. Hence, the TETRAD_TYPE is 
set to zero, the base vectors BASE_0 to BASE_3 build the identity matrix, and there is no boost 
transformation (the last value is the boost velocity). Since we use geometrical units, the speed 
of light and the gravitational constant are set to unity. 
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Appendix B. Periodic orbits by Levin and Perez-Giz 

The periodic orbits in the Schwarzschild space-time shown in the paper by Levin and Perez- 
Giz |9 1, can be reproduced in isotropic coordinates in the following way. For a nearly arbitrary 
radial coordinate r, the initial velocity /3 and the initial direction ^ read 



^ arcsin^ ^ (B.l) 

The isotropic position is given by 

X = ^ (2r - r, + 2 ^Jr{r~r,)) . (B.2) 

For example, to reproduce the periodic orbit of Fig. 12 in their paper with parameters L — 3>.9 
and E = 0.987160, we set r = 10 and obtain j3 w 0.423147, t, « 56.62467, andx « 8.972136. 
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